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Abstract 

In this paper we combine finite difference approximations (for spatial derivatives) and 
collocation techniques (for the time component) to numerically solve the two dimensional 
heat equation. We employ respectively a second-order and a fourth-order schemes for the 
spatial derivatives and the discretization method gives rise to a linear system of equations. 
We show that the matrix of the system is non-singular. Numerical experiments carried 
out on serial computers, show the unconditional stability of the proposed method and 
the high accuracy achieved by the fourth-order scheme. 
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1 Introduction 


We consider the two dimensional heat equation: 

du , x 2 (d 2u , d 2 u f A 

m (x ' y ' t) = “ V x,y ’ 7 ’ 

' t ’ 

where fi = [0, 1] x [0, 1], and with the initial condition 


{x,y,t) G 0 x [0,oo) (1) 


u(x,y,0) -ip{x,y), (x,y)<=Cl, 


and the boundary conditions 

u(0,j/,t) = fo{y,t), = /i(y,i), u(x,0,t) = g 0 (x,t), and «(x, l,t) = gi{x,t) f or t > 0. 
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We assume that /o, /i, go and g\ are smooth functions in the variable i, i.e., their first 
derivatives with respect to t exist and axe continuous. 

If (1) is discretized with standard or high-order finite difference approximations, the 
resulting method leads to a stability condition. For a reasonable number of mesh points in 
the spatial direction, this typically requires a very small time step to satisfy the stability 
requirement. Even if the approximations produce implicit methods, the computational com- 
plexity considerably increases, especially if high-order formulas are employed [5]. In addition, 
these techniques, when implemented on parallel computers can only allow the parallelization 
in space, i.e., at each time step, spatial grid points are partitioned and assigned to processors; 
the solution is then computed before we move to the next time step. 

Recently, Jezequel [3] combined the standard finite difference approximation for the 
spatial derivative and collocation technique for the time component to numerically solve the 
one dimensional heat equation. The method (called implicit collocation method) is uncon- 
ditionally stable. Its principle is as follows: after discretization in space of the problem, the 
solution is approximated at each spatial grid point by a polynomial depending on time. The 
resulting derivation produces a linear system of equations. The order of the method is in 
space the order of difference approximation and in time the degree of the polynomial. 

In this paper, we extend Jezequel’s work [3] to the two dimensional heat equation. In 
addition to the spatial discretization with the standard second-order formula, we also present 

T ’ 

discretization based on a fourth-order formula. For the two formulas, we show that the 
matrix arising from the system of equations is non-singular and we present their respective 
accuracies. Our numerical experiments are carried out on a serial computer. 

An outline of the paper is as follows. In Section 2, we explain the basic principle behind 
the implicit collocation method. Section 3 presents the derivation of the system of linear 
equations when the fourth-order and the second-order formulas are respectively utilized for 
the spatial derivatives. Numerical results axe given in Section 4. In Section 5, we discuss some 
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issues related to the implicit collocation method. Finally, some conclusions are formulated in 
Section 6. 

2 Principle of the Implicit Collocation Technique 

The idea behind the technique can be described as follows: 

1. We start with a time dependent partial differential equation (PDE). 

2. The PDE is discretized in space, giving rise to a system of ordinary differential equations 
with unknown functions at each spatial grid point. 

3. The implicit collocation method consists of approximating at each spatial grid point 
the solution by a polynomial that depends on time. To solve the PDE with implicit 
collocation method is to determine the coefficients of all polynomials. 

4. Depending of the PDE, we obtain a linear or non-linear system of equations (where the 
unknowns are the coefficients) that can be solved by a direct or iterative method. 

5. Once the coefficients of the polynomials are determined, the -approximated solution 

of the PDE is computed on a given time interval that depends on the degree of the 
polynomials. , * 

One of the main advantages of the implicit collocation method is that if it is efficiently 
implemented on distributed memory computers, the parallelization is carried out both across 
time and space [3]. 

In the next section we use this description to derive the system of equations. The 
main presentation focuses on the fourth-order finite difference approximations for spatial 
derivatives. 
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3 Derivation of the System of Equations 

3.1 Discretization Procedure 

Let h = l / n and At be the spatial and time mesh-widths respectively. We can subdivide the 
spatial domain and consider the time step as follows: 

Xj = ih, yj — jh) j 0, 1, . . . ,n 

tk — kAt , k = 0 , 1 , . . « 

For simplicity we write the approximated solution of u and its time derivative at the spatial 
grid points {xi>yj) as: 

Ui ,j{t) = u{xi,yj,t), and U-j(t) = — {x u yj,t). 

At any given time t , if we use the discretization of the steady state Poisson equation 
with a fourth-order (FO) scheme [1], we can approximate the spatial derivatives of (1). We 
obtain for any grid point ( Xi,yj ) , *, j = 1,. . . ,n — 1: 

\ KuM + By+1 M + UU,i(t) + Olj-M + 8^(0] 

= ^(Vi+uW+Oirf+iW + Oi-ijW + Oy-iW) (2) 

+ (Ph-ij+iM + Ui-\j+\(t) -I- + Pj+ij_i(<)) — 20Uij(t)], 

with the conditions > , 

U%,j(0) = ip{x i) y j ), j '■ 

Uoj{t) = fo{yj,t), U 0>j (t) = -^-(yj,t), 

Unj{t) = MVj , <). C^(0 = *), 

Uifi(t) = ^ii.oCO = ^ 

= <7l (*»,*), 

Eq. 2 is a system of (n - 1) x (n - 1) ordinary differential equations and for any value 
of t, it is fourth-order in space. 
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Remark 1 In case the spatial derivative in (1) is discretized with the standard second-order 
(SO) finite difference approximation , we obtain the system: 

Ulj(t) = ^ \U i+ ij{t) + Ui, j+l (t) + Ui-uit) + Uij-iit) - 4U t< M . (3) 

Here the conditions on the time derivatives on boundary points are not employed . 


Now it remains to introduce the concept of implicit collocation methods in our deriva- 
tions. 

Let Pi,j{t) be the polynomial of degree r satisfying the system (2) at the spatial grid 
point (xi,yj) and at times t^ ~ kAt (fc = 0, . . . ,r — 1). Then for any i,j = 1, ... ,n — 1 and 
k = 0, . . . , r — 1 , we have 

— a hj,T^k "I" a tj,r- ll'k + * * * + Qijyltk + a iJ, 0- 

The coefficients dij t o are determined from the initial condition: 

diyjyQ ~ Pi j(0) = = 

To solve the system (2) by the collocation method is to determine the coefficients a t • • • > a ij y r> 
for i, j = 1, . . . ,n — 1. After some algebraic manipulations (see [3] for details on the one di- 
mensional heat equation) we obtain the linear system of r x (n — 1) x (n — 1) equations 

AX = S, (4) 

) 

where A is a block-tridiagonal matrix given by ( 

i ' 

A = tri [Aj-i, Aj, • 


Aj- 1 , Aj and Aj+ 1 are square matrices (with r x (n — 1) rows) defined as 

1 h 2 

Aj.r = tri -E,-^E'-4E,~ 

"l h2 1 /i2 

Aj = tri -^E' -4E,4~E' + 20E,-— E'-4E 

J 2 or a z 2 or . 

J n— 1 

1 h 2 

A j+X = tri -E y -—^E'-4E-E . 

J n-l 
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The subscript n — 1 determines the number of block-rows- E and E f are r x r nonsymmetric 


matrices 


E = 

( *5 C 1 • 

*1 <i _1 ' 

k \ 
• h 

, E' = 

(rf 0 ~' (r-lW 2 ' 

rt\ 1 (r - 1)*; 2 • 

• 2f 0 1 \ 

• 2ii 1 


\ tr-\ tr-l ' 

■ ■ tr—l / 


Uci (r-l)C? ■' 

’ ■ 2< r _i 1 J 


the vector X of r x (n — 1) x (n ~ 1) unknowns and the right hand side S are 


( ( °U,r \ \ 


X = 



\ &n— 2,n— 1,1 
a n— l,n— l,r \ 


\ V <3n-l,n-l,r / / 


( 


s = 


( 4(02,1,0 + fl i,2,o) + ^2,2,0 — 20ai ? i,o + 4(?7i ) o(^o) + £^o,i(^o)) j 

+Uo,2(to) + Uo,o(to) + U 2 fi (to) - + U^(to)) | 

4(a2,l,0 + ® l,2,o) + 02,2,0 — 2001,1,0 + 4({7i, 0 (t r _i + f^0,l(^r-l) 

^ +t^0,2(<0) + Co.o(^r-l) + ^2,o(*r-l) ~ 2 £*(^l,o( 4 r-l) + Co,i(*r-l)) ) 

/ 4(oj + ij,o 4- Ojj+ 1,0 + Oj_i,j,o +r a t,i-i,o) + a «+ij+l,0 \ 

+Oi-ij+i,o + Oj_i,j_i,o + a,+i,j-i,o — 20ojj,o 


I 4(aj + ij,o + Oij+ 1,0 + °t-ij,o + atj- 1 , 0 ) + a «+iJ+i,o ! 
\ +o*-ij+i,o + + a»+ij-i,o - 20ojj,o / 


( 4(a n _2,n— 1,0 + O n _i,n_2,o) + a n- 2,n-2,0 — 20a„_i, tv _i,o + 4(£/n_i, n (*o) + C n , n . 

+£ 7 n,n-2(*o) + U ntn {to) + U n -2, n (to) ~ 2 (^»-l,n(*o) + ^n.n-l^o)) 

4(a n _2,n-l,0 + 011 - 1 , 71 - 2 , 0 ) + O n _2,n— 2,0 — 20a n _i, n _i,o + 4(?/ n _l, n (*r-l) + U n>n . 

^ +C^n,n— 2(^r—l) 4" Cn,n(^r— l) + Un-2,n{tr-l) — 2 ^(^n-l,n(^r-l) + E nn _i {t 


-i(*o)) \ 

-l(<r- 1)) 

-l)) / / 
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Remark 2 With the SO scheme, the matrix A is still block-tridiagonal and 


Aj-\ 

Aj 

Aj + 1 


diag\-E ]„_i, 


= tri 


-E,%E' + 4E,-E 
or 


J n— 1 


= diag[-E\ n -i. 


The vector of unknowns remains the same and the right hand side S does not involve time 
derivatives at boundary points . 

Remark 3 A is a matrix with bandwith equal to (2n — 2)r and (2 n+ l)r for the SO and FO 
spatial schemes respectively. The block structure of the matrix A for the SO or FO scheme , is 
similar to the one obtained from the discretization of the two dimensional steady-state Poisson 
equation with the SO or FO scheme . In the latter , instead of having the block matrices E and 
E f , we have constant coefficients. 

Remark 4 To obtain the solution , the coefficients (i, j = 1, . . ■ ,n — 1 and k = 1, . . . } r) 
are first evaluated and then the approximated values (hj — — 1 and k = 

1 } . . . ,r - l) are calculated . Each of these two steps can be carried out in parallel In addi- 
tion, the collocation method does not only consist of determining the only , but also 

calculates coefficients of polynomials approaching the solution in the interval [to^r-i]- 

/ ■ 

3.2 Non- Singularity of the Matrix A 

If we focus on the FO spatial approximation, Eq. 4 lias a solution only if the matrix A is 
non-singular. To determine the non-singularity of ^4, it is enough to study its block-diagonal 
elements. 

Let D = E f + 20E be the block-diagonal of the matrix A{ defined above. We can 
state the following theorem. 

Theorem 1 The determinant of the matrix D is zero if and only if two values t[- \ and i m -\ 
are identical 
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Proof: D is a r x r matrix where 


Did = + 1 -i) + 20 *i— lj t r i- 1 > 0, ij = 1 

For a given i, it is impossible that Djj = 0 for all values of j. Indeed the cofficients of the 
last column of the matrix D are all non zero: D Jir = + 20^_i ^ 0. D can be seen as a 

Vandermonde- like matrix. Since all the 1 %- 1 are all different, D is a non-singular matrix, n 

Conjecture 1 The determinant of the matrix A is non zero and the linear system (4) has 
a unique solution . 

The same conclusion can be obtained if we employ the SO spatial scheme. 

4 Numerical Experiments 

Consider Eq. 1 with the conditions 

a = ti(x,y,0) = sin7rir + sin7ry, 

7 r 

fo(y,t) = fi(y,t) = e -t sinny, go(x,t) = gi(x,t) = e -t sin7rx. 

The exact solution is given by u(x } y,t) = (sinTrx + sin7ry)e~ < . 

For the second-order (SO) and the fourth-order (FO) spatial discretizations respectively, 
the implicit collocation technique was implemented an SGI 02 Workstation in Fortran 77. 
To solve the linear system of equations, we used the decojnposition algorithm for inverting 
asymmetric and indefinite matrices proposed by Luo [4]. The method is easy to programmed, 
requires only the storage of the matrix and the right hand side of Eq. 4 and was chosen for 
its parallelization potential. 

To test the accuracy and stability of the implicit collocation technique, we present for 
= 0.1 and At = 0.01 and for different values of r and n, the absolute maximum error 
(obtained by comparing the true solution with the approximated one) achieved in the interval 
[0, (r — 1)A<]. We compare errors obtained with the second-order (SO) and the fourth-order 
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(FO) spatial schemes. We do not report the timing results because the two schemes give the 
same elapsed time. This is due to the fact that in the implementation of Luo’s algorithm, we 
did not take advantage of the respective matrix bandwidths; even if we did, the difference in 
elapsed times would be negligible since their bandwidths are comparable. 

To provide an idea on how large the system of equations is, we first present in Table 1, 
the dimension of the matrix A in (4) for different values of n and r. 


n 

r = 3 

r = 4 

r = 5 

4 

8 

16 

27 x 27 
147 x 147 
675 x 675 

36 x 36 
196 x 196 
900 x 900 

45 x 45 
245 x 245 
1125 x 1125 


Table 1: Dimension of the matrix A for different values of n and r. 

We report in Tables 2, 3 and 4 absolute maximum errors for r = 3,4,5 respectively 
when the spatial mesh-width varies. The maximum error in the interval [0, (r — 1) At] was 



At = 

= 0.1 

At = 

0.01 

n 

SO 

FO 

SO 

FO 

4 

1.61 (-2) 

5.35(-4) 

1.97(-3) 

6.27(-5) 

8 

4.14(-3) 

3.33(-5) 

5.01 (-4) 

3.90(-6) 

16 

1.04(-3) 

1.89(-6) 

1.25 (-4) 

2.43(-7) 


Table 2: Maximum error obtained with the second-order (SO) and fourth-order (FO) spatial 
schemes for different values of At when r = 3. 

obtained for t = (r — 1) At and the error increases as t does within the interval. This result is 
consistent with the one achieved by Jezequel while solving the one dimensional heat equation 
using the second order scheme and collocation technique [3]. 

We observe that for given r and At, SO indeed produces solutions of second-order 
accuracy and FO of fourth-order accuracy. In addition, with FO the accuracy is far more 
better. These findings are consistent with the ones obtained by Gupta et al [2] and Zhang 
et al [6] while solving the steady state Poisson equation with the two schemes. 
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At = 

= 0.1 

At = 


n 

SO 

FO 

SO 

FO 

4 

2.14(-2) 

7.10(-4) 

2.93(-3) 

9.35(-5) 

8 

5.52(-3) 

4.28(-5) 

7.44(-4) 

5.80(-6) 

16 

1.39(-3) 

2.16(-6) 

1.86(-4) 

3.61(-7) 


Table 3: Maximum error obtained with the second-order (SO) and fourth-order (FO) spatial 
schemes for different values of At when r — 4 . 



At = 

= 0.1 

At = 

0.01 


SO 

FO 

SO 

FO 


2.51 (-2) 

8.34(-4) 

3.86(-3) 

1.23(-4) 

8 

6.49(-3) 

5.10(-5) 

9.82(-4) 

7.66(-6) 

16 

1.63(-3) 

3.18(-6) 

2.46(-4) 

4.76(-7) 


Table 4: Maximum error obtained with the second-order (SO) and fourth-order (FO) spatial 
schemes for different values of At when r = 5. 

5 Discussion 

5.1 Degree of the Polynomials 

The order of the implicit collocation technique is the order of the difference scheme in space 
and the order of the polynomial in time (namely r). The question that arises is how do we 
choose r in order to obtain high accurate approximated solutions of (1)? By increasing the 
value of r, we dot not only increase the lenght of the time interval where the solution is to 
be found but also the size of the linear system of equations. In [3], Jezequel studied the 
numerical validity of the coefficients of the polynomials. She defined the optimal degree r op t 
of the polynomials to be the highest integer for which all the coefficients of the polynomials 
and the approximated solution remain significant. She found that the degree r op t increases 
as the time and spatial mesh sizes increase but r op t is not arbitrary large. 
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5.2 Merit of the Method 


Here we attempt to computationally compare the implicit collocation method and other 
standard implicit methods (arising from classical finite difference approximations in both 
space and time). 

Assume that the spatial domain has (n — 1) x (n — 1) interior grid points. The derivations 
with a standard implicit method (SI) produces a system of equations with N S { = (n — 1) x 
(n — 1) unknowns whereas with the implicit collocation (IC) method, we obtain a system 
with Ni C = r x (n - 1) x (n - 1) unknowns, where r is the degree of the polynomials. 

Let A ti C be the time step used for IC. The solution using IC, can be determined at 
any point in the interval [0, (r — 1) A i$ c ]. Let (r — l)m be the number of equidistant points 
where the solution is to be computed in this time interval. To determine the solution at 
the same points of the interval with SI, (r — l)m time iterations must be carried out. The 
corresponding time step is A t s i = A U c jm. 

The implementation of IC requires the solution of a linear system of Nj c equations; its 
computational cost is Ci c ~ Nf c . For SI, we need to find the solution of a linear system of 
N S { equations (r — l)m times; the cost is then C s i « (r — l)mAT^. If we assume that Ci c and 
C S i are equal, then r 3 = (r — l)m orm = r 2 + r + l + l/(r-l). We can conclude that if 
m <r 2 + r + 2, then SI is computationally less expensive than IC, and ifm>r 2 -f-r + 2,IC 
is cheaper. 

To summarize, the implicit collocation method is benefitial with respect to standard 
implicit methods, if for given r and A t, the number of equidistant time points in the interval 
[0, (r — 1) At] is at least greater than r 3 + r. 

6 Conclusions 

We have carried out the numerical approximation of the two dimensional heat equation by 
using finite difference schemes (for spatial derivatives) and implicit collocation technique (for 
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time component). The proposed method is unconditionally stable. We employed second- 
order and a fourth-order spatial schemes respectively and showed that both have the same 
computational complexity and that the fourth-order one clearly produces more accurate so- 
lutions. 

The main advantage of the implicit collocation technique is not only its stability condition 
but also the fact that it can be implemented on distributed memory computers where the 
parallelization strategy is performed both across time and space. In future works, we plan 
to implement the method on parallel computers and to extend our analysis to the three 
dimensional heat equation. 
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